package scu.maqiang.moleculardynamics;

import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.MVO;

import java.io.*;
import java.security.SecureRandom;
import java.util.*;

/**
 * 分子动力学类
 */
public class MolecularDynamics {
    public MolecularDynamics(int nDim) {
        this.nDim = nDim;
    }

    public void initializePosition(String fileName) {
        try (BufferedReader br = new BufferedReader(new FileReader(fileName))) {
            String ll = br.readLine();
            number = Integer.parseInt(ll);
            String[] data = br.readLine().split("\\s+");
            boxMatrix = MVO.zero(nDim, nDim);
            if (nDim == 3) {
                if (data.length == 3) {
                    boxMatrix[0][0] = Double.parseDouble(data[0]);
                    boxMatrix[1][1] = Double.parseDouble(data[1]);
                    boxMatrix[2][2] = Double.parseDouble(data[2]);
                } else if (data.length == 9) {
                    boxMatrix[0][0] = Double.parseDouble(data[0]);
                    boxMatrix[0][1] = Double.parseDouble(data[1]);
                    boxMatrix[0][2] = Double.parseDouble(data[2]);
                    boxMatrix[1][0] = Double.parseDouble(data[3]);
                    boxMatrix[1][1] = Double.parseDouble(data[4]);
                    boxMatrix[1][2] = Double.parseDouble(data[5]);
                    boxMatrix[2][0] = Double.parseDouble(data[6]);
                    boxMatrix[2][1] = Double.parseDouble(data[7]);
                    boxMatrix[2][2] = Double.parseDouble(data[8]);
                } else {
                    System.out.println("The second line should have 3 or 9 items in 3-dimension.");
                    System.exit(1);
                }
                boxMatrixInv = MVO.inverseMat3(boxMatrix);
            } else if (nDim == 2) {
                if (data.length == 2) {
                    boxMatrix[0][0] = Double.parseDouble(data[0]);
                    boxMatrix[1][1] = Double.parseDouble(data[1]);
                } else if (data.length == 4) {
                    boxMatrix[0][0] = Double.parseDouble(data[0]);
                    boxMatrix[0][1] = Double.parseDouble(data[1]);
                    boxMatrix[1][0] = Double.parseDouble(data[2]);
                    boxMatrix[1][1] = Double.parseDouble(data[3]);
                } else {
                    System.out.println("The second line should have 2 or 4 items in 2-dimension.");
                    System.exit(1);
                }
                boxMatrixInv = MVO.inverseMat2(boxMatrix);
            }
            thickness = MVO.getThickness(boxMatrix);
//            System.out.println(Arrays.toString(thickness));
//            System.exit(1);

            NN = new int[number];
            NL = new int[number * MN];
            mass = new double[number];
            x0 = new double[number];
            y0 = new double[number];
            z0 = new double[number];
            x = new double[number];
            y = new double[number];
            z = new double[number];
            vx = new double[number];
            vy = new double[number];
            vz = new double[number];
            fx = new double[number];
            fy = new double[number];
            fz = new double[number];
            type = new String[number];

            System.out.println("box matrix = ");
            for (int i = 0; i < 3; ++i) {
                for(int j = 0; j < 3; ++j) {
                    System.out.print(boxMatrix[i][j] + "  ");
                }
                System.out.println();
            }

            System.out.println("box matrix inverse = ");
            for (int i = 0; i < 3; ++i) {
                for(int j = 0; j < 3; ++j) {
                    System.out.print(boxMatrixInv[i][j] + "  ");
                }
                System.out.println();
            }

            for(int i = 0; i < number; i++) {
                String[] tokens =  br.readLine().split("\\s+");
                type[i] = tokens[0];
                x[i] = Double.parseDouble(tokens[1]);
                y[i] = Double.parseDouble(tokens[2]);
                z[i] = Double.parseDouble(tokens[3]);
                mass[i] = Double.parseDouble((tokens[4]));
            }
        } catch (IOException e) {
            e.printStackTrace();
        }

    }

    public void initializeTimeIntegrationInfo(String fileName) {

        try (BufferedReader br = new BufferedReader(new FileReader(fileName))) {
            String line;
            while((line = br.readLine()) != null) {
                String[] tokens = line.split("\\s+");
                if (tokens.length > 0) {
                    if (Objects.equals(tokens[0], "time_step")) {
                        timeStep = Double.parseDouble(tokens[1]);
                        if (timeStep < 0) {
                            System.out.println("timeStep should >= 0.");
                            System.exit(1);
                        }
                        System.out.println("timeStep: " + timeStep);
                    } else if (Objects.equals(tokens[0], "run")) {
                        numSteps = Integer.parseInt(tokens[1]);
                        if (numSteps < 1) {
                            System.out.println("numSteps should >= 1.");
                            System.exit(1);
                        }
                        System.out.println("numStep: " + numSteps);
                    } else if (Objects.equals(tokens[0], "velocity")) {
                        temperature = Double.parseDouble(tokens[1]);
                        if (temperature < 0) {
                            System.out.println("temperature >= 0.");
                            System.exit(1);
                        }
                        System.out.println("temperature: " + temperature);
                    } else if (Objects.equals(tokens[0], "neighbor_flag")) {
                        neighbor_flag = Integer.parseInt(tokens[1]);
                        if (neighbor_flag < 0 || neighbor_flag > 2) {
                            System.out.println("neighbor_flag can only be 0 or 1 or 2.");
                            System.exit(1);
                        }
                        System.out.println("neighbor_flag: " + neighbor_flag);
                    } else if (tokens[0].charAt(0) != '#') {
                        System.out.println(tokens[0] + " is not a valid keyword.");
                        System.exit(1);
                    }
                }
            }
            timeStep /= TIME_UNIT_CONVERSION;
        } catch (Exception e) {
            throw new RuntimeException(e);
        }
    }

    public void initializeUniformVelocity(double a, double b) {
        Random rand = new Random(System.nanoTime());
        //质心速度
        double[] centerOfMassVelocity = {0.0, 0.0, 0.0};
        double totalMass = 0.0;
        for (int n = 0; n < number; ++n) {
            totalMass += mass[n];
            vx[n] = a + rand.nextDouble() * (b - a);
            vy[n] = a + rand.nextDouble() * (b - a);
            vz[n] = a + rand.nextDouble() * (b - a);
            centerOfMassVelocity[0] += mass[n] * vx[n];
            centerOfMassVelocity[1] += mass[n] * vy[n];
            centerOfMassVelocity[2] += mass[n] * vz[n];
        }
        centerOfMassVelocity[0] /= totalMass;
        centerOfMassVelocity[1] /= totalMass;
        centerOfMassVelocity[2] /= totalMass;
        for (int n = 0; n < number; ++n) {
            vx[n] -= centerOfMassVelocity[0];
            vy[n] -= centerOfMassVelocity[1];
            vz[n] -= centerOfMassVelocity[2];
        }
        scaleVelocity();
    }

    protected void scaleVelocity() {
        final double computedTemperature = computeKineticEnergy() * 2.0 / (3.0 * K_B * number);
        double scaleFactor = Math.sqrt(temperature / computedTemperature);
        for (int n = 0; n < number; ++n) {
            vx[n] *= scaleFactor;
            vy[n] *= scaleFactor;
            vz[n] *= scaleFactor;
        }
    }

    protected double computeKineticEnergy() {
        double kineticEnergy = 0.0;
        for (int n = 0; n < number; ++n) {
            kineticEnergy += mass[n] * (vx[n] * vx[n] + vy[n] * vy[n] + vz[n] * vz[n]);
        }
        return kineticEnergy * 0.5;
    }

    public double applyMicOne(double x12) {
        double result = x12;
        if (x12 < -0.5) {
            result =  x12 + 1.0;
        } else if (x12 > +0.5) {
            result = x12 - 1.0;
        }
        return result;
    }

    public double[] applyMic(double x12, double y12, double z12) {
        double sx12 = boxMatrixInv[0][0] * x12 + boxMatrixInv[0][1] * y12 + boxMatrixInv[0][2] * z12;
        double sy12 = boxMatrixInv[1][0] * x12 + boxMatrixInv[1][1] * y12 + boxMatrixInv[1][2] * z12;
        double sz12 = boxMatrixInv[2][0] * x12 + boxMatrixInv[2][1] * y12 + boxMatrixInv[2][2] * z12;
        sx12 = applyMicOne(sx12);
        sy12 = applyMicOne(sy12);
        sz12 = applyMicOne(sz12);
        double[] newXYZ = {0.0, 0.0, 0.0};
        newXYZ[0] = boxMatrix[0][0] * sx12 + boxMatrix[0][1] * sy12 + boxMatrix[0][2] * sz12;
        newXYZ[1] = boxMatrix[1][0] * sx12 + boxMatrix[1][1] * sy12 + boxMatrix[1][2] * sz12;
        newXYZ[2] = boxMatrix[2][0] * sx12 + boxMatrix[2][1] * sy12 + boxMatrix[2][2] * sz12;
        return newXYZ;
    }

    public double applyPbcOne(double x) {
        double result = x;
        if (x < 0.0) {
            result =  x + 1.0;
        } else if (x > 1.0) {
            result = x- 1.0;
        }
        return result;
    }

    public void applyPbc() {
        for(int n = 0; n < number; ++n) {
            double sx = boxMatrixInv[0][0] * x[n] + boxMatrixInv[0][1] * y[n] + boxMatrixInv[0][2] * z[n];
            double sy = boxMatrixInv[1][0] * x[n] + boxMatrixInv[1][1] * y[n] + boxMatrixInv[1][2] * z[n];
            double sz = boxMatrixInv[2][0] * x[n] + boxMatrixInv[2][1] * y[n] + boxMatrixInv[2][2] * z[n];
            sx = applyPbcOne(sx);
            sy = applyPbcOne(sy);
            sz = applyPbcOne(sz);
            x[n] = boxMatrix[0][0] * sx + boxMatrix[0][1] * sy + boxMatrix[0][2] * sz;
            y[n] = boxMatrix[1][0] * sx + boxMatrix[1][1] * sy + boxMatrix[1][2] * sz;
            z[n] = boxMatrix[2][0] * sx + boxMatrix[2][1] * sy + boxMatrix[2][2] * sz;
        }
    }

    public void computeLJForce() {
	    final double epsilon = 1.032e-2;
        final double sigma = 3.405;
        final double cutoff = 9.0;
        final double cutoffSquare = cutoff * cutoff;
        final double sigma3 = sigma * sigma * sigma;
        final double sigma6 = sigma3 * sigma3;
        final double sigma12 = sigma6 * sigma6;
        final double e24s6 = 24.0 * epsilon * sigma6;
        final double e48s12 = 48.0 * epsilon * sigma12;
        final double e4s6 = 4.0 * epsilon * sigma6;
        final double e4s12 = 4.0 * epsilon * sigma12;
        pe = 0.0;
        Arrays.fill(fx, 0.0);
        Arrays.fill(fy, 0.0);
        Arrays.fill(fz, 0.0);
        if (neighbor_flag == 0) {
            for (int i = 0; i < number - 1; ++i) {
                for (int j = i + 1; j < number; ++j) {
                    double xij = x[j] - x[i];
                    double yij = y[j] - y[i];
                    double zij = z[j] - z[i];
                    double[] newXYZ = applyMic(xij, yij, zij);
                    xij = newXYZ[0];
                    yij = newXYZ[1];
                    zij = newXYZ[2];
                    final double r2 = xij * xij + yij * yij + zij * zij;
                    if (r2 > cutoffSquare)
                        continue;

                    final double r2inv = 1.0 / r2;
                    final double r4inv = r2inv * r2inv;
                    final double r6inv = r2inv * r4inv;
                    final double r8inv = r4inv * r4inv;
                    final double r12inv = r4inv * r8inv;
                    final double r14inv = r6inv * r8inv;
                    final double f_ij = e24s6 * r8inv - e48s12 * r14inv;
                    pe += e4s12 * r12inv - e4s6 * r6inv;
                    fx[i] += f_ij * xij;
                    fx[j] -= f_ij * xij;
                    fy[i] += f_ij * yij;
                    fy[j] -= f_ij * yij;
                    fz[i] += f_ij * zij;
                    fz[j] -= f_ij * zij;
                }
            }
        } else {
            for(int i = 0; i < number; i++) {
                final double xi = x[i];
                final double yi = y[i];
                final double zi = z[i];
                for (int jj = 0; jj < NN[i]; ++jj) {
                    final int j = NL[i * MN + jj];
                    double xij = x[j] - xi;
                    double yij = y[j] - yi;
                    double zij = z[j] - zi;
                    double[] newXYZ = applyMic(xij, yij, zij);
                    xij = newXYZ[0];
                    yij = newXYZ[1];
                    zij = newXYZ[2];
                    final double r2 = xij * xij + yij * yij + zij * zij;
                    if (r2 > cutoffSquare)
                        continue;

                    final double r2inv = 1.0 / r2;
                    final double r4inv = r2inv * r2inv;
                    final double r6inv = r2inv * r4inv;
                    final double r8inv = r4inv * r4inv;
                    final double r12inv = r4inv * r8inv;
                    final double r14inv = r6inv * r8inv;
                    final double f_ij = e24s6 * r8inv - e48s12 * r14inv;
                    pe += e4s12 * r12inv - e4s6 * r6inv;
                    fx[i] += f_ij * xij;
                    fx[j] -= f_ij * xij;
                    fy[i] += f_ij * yij;
                    fy[j] -= f_ij * yij;
                    fz[i] += f_ij * zij;
                    fz[j] -= f_ij * zij;
                }
            }
        }
    }

    public void integrate(final boolean isStepOne, final double timeStep) {
        //计算半个时间步长
	    final double timeStepHalf = timeStep * 0.5;
        for (int n = 0; n < number; ++n) {
            final double mass_inv = 1.0 / mass[n];
            final double ax = fx[n] * mass_inv;
            final double ay = fy[n] * mass_inv;
            final double az = fz[n] * mass_inv;
            vx[n] += ax * timeStepHalf;
            vy[n] += ay * timeStepHalf;
            vz[n] += az * timeStepHalf;
            if (isStepOne) {
                x[n] += vx[n] * timeStep;
                y[n] += vy[n] * timeStep;
                z[n] += vz[n] * timeStep;
            }
        }
    }


    public void findNeighbor_ON2() {
        final double cutoffSquare = cutoffNeighbor * cutoffNeighbor;
        Arrays.fill(NN, 0);
        for (int i = 0; i < number - 1; ++i) {
            final double x1 = x[i];
            final double y1 = y[i];
            final double z1 = z[i];
            for (int j = i + 1; j < number; ++j) {
                double xij = x[j] - x1;
                double yij = y[j] - y1;
                double zij = z[j] - z1;
                double[] newXYZ = applyMic(xij, yij, zij);
                xij = newXYZ[0];
                yij = newXYZ[1];
                zij = newXYZ[2];
                final double distanceSquare = xij * xij + yij * yij + zij * zij;
                if (distanceSquare < cutoffSquare) {
                    NL[i * MN + NN[i]++] = j;
                    if (NN[i] > MN) {
                        System.out.println("Error: number of neighbors for atom " + i + " exceeds " + MN);;
                        System.exit(1);
                    }
                }
            }
        }
    }

    public void findNeighbor_ON2_2() {
        neighbor_list = new ArrayList[number];
        for(int i = 0; i < number; i++) {
            neighbor_list[i] = new ArrayList<Integer>();
        }
        final double cutoffSquare = cutoffNeighbor * cutoffNeighbor;
        for (int i = 0; i < number - 1; ++i) {
            final double x1 = x[i];
            final double y1 = y[i];
            final double z1 = z[i];
            for (int j = i + 1; j < number; ++j) {
                double xij = x[j] - x1;
                double yij = y[j] - y1;
                double zij = z[j] - z1;
                double[] newXYZ = applyMic(xij, yij, zij);
                xij = newXYZ[0];
                yij = newXYZ[1];
                zij = newXYZ[2];
                final double distanceSquare = xij * xij + yij * yij + zij * zij;
                if (distanceSquare < cutoffSquare) {
                    neighbor_list[i].add(j);
                }
            }
        }
    }

    public void updateCoords() {
        for (int n = 0; n < number; ++n) {
            x0[n] = x[n];
            y0[n] = y[n];
            z0[n] = z[n];
        }
    }

    public boolean checkIfNeedUpdateNeighbor() {
        boolean needUpdate = false;
        for (int n = 0; n < number; ++n) {
            double dx = x[n] - x0[n];
            double dy = y[n] - y0[n];
            double dz = z[n] - z0[n];
            if (dx * dx + dy * dy + dz * dz > 0.25) {
                needUpdate = true;
                break;
            }
        }
        return needUpdate;
    }

    public void findCell(final double[] r, double cutoffInverse, int[] numCells, int[] cell) {
        //首先将其变换到标准单位正方体坐标中
        double[] s = new double[3];
        s[0] = boxMatrixInv[0][0] * r[0] + boxMatrixInv[0][1] * r[1] + boxMatrixInv[0][2] * r[2];
        s[1] = boxMatrixInv[1][0] * r[0] + boxMatrixInv[1][1] * r[1] + boxMatrixInv[1][2] * r[2];
        s[2] = boxMatrixInv[2][0] * r[0] + boxMatrixInv[2][1] * r[1] + boxMatrixInv[2][2] * r[2];
        //找到在各个维度上的索引
        for (int d = 0; d < 3; ++d) {
            cell[d] = (int) Math.floor(s[d] * thickness[d] * cutoffInverse);
            //如果小于0, 则放到前一个索引中
            if (cell[d] < 0) {
                cell[d] += numCells[d];
            }
            //如果大于等于最大的索引值, 则放到最后一个索引中
            if (cell[d] >= numCells[d]) {
                cell[d] -= numCells[d];
            }
        }
        //计算得到整体索引值
        cell[3] = cell[0] + numCells[0] * (cell[1] + numCells[1] * cell[2]);
        //System.out.println(Arrays.toString(cell));
    }


    public void findNeighbor() {
        if (checkIfNeedUpdateNeighbor()) {
            numUpdates++;
            applyPbc();
            if (neighbor_flag == 1) {
                findNeighbor_ON1();
            } else if (neighbor_flag == 2) {
                findNeighbor_ON2();
            }
            updateCoords();
        }
    }

    public void findNeighbor_ON1() {
        final double cutoffInverse = 1.0 / cutoffNeighbor;
        double cutoffSquare = cutoffNeighbor * cutoffNeighbor;

        int[] numCells = new int[4];
        for (int d = 0; d < 3; ++d) {
            numCells[d] = (int) Math.floor(thickness[d] * cutoffInverse);
        }
        numCells[3] = numCells[0] * numCells[1] * numCells[2];
        // cout << "numCells" << endl;
        // cout << numCells[0] << "  " << numCells[1] << "  " << numCells[2] << endl;
        // exit(1);
        int[] cell = new int[4];

        int[] cellCount = new int[numCells[3]];

        int[] cellCountSum = new int[numCells[3]];

        //cout << "OKOKOK" << endl;
        //cout << "Thickness" << endl;
        //cout << thickness[0] << endl;
        //cout << thickness[1] << endl;
        //cout << thickness[2] << endl;

        //统计这些每个Cell中有多少个原子
        for (int n = 0; n < number; ++n) {
            final double[] r = {x[n], y[n], z[n]};
            findCell(r, cutoffInverse, numCells, cell);
            ++cellCount[cell[3]];
        }
        // cout << "cellCount: " << endl;
        // for(int i = 0; i < numCells[3]; ++i) {
        //     std::cout << cellCount[i] << endl;
        // }
        // cout << "OKOKOK" << endl;

        for (int i = 1; i < numCells[3]; ++i) {
            cellCountSum[i] = cellCountSum[i - 1] + cellCount[i - 1];
        }

        Arrays.fill(cellCount, 0);

        int[] cellContents = new int[number];
        //cout << "OKOKOK" << endl;
        for (int n = 0; n < number; ++n) {
            final double[] r = {x[n], y[n], z[n]};
            findCell(r, cutoffInverse, numCells, cell);
            //cout << "--" << endl;
            cellContents[cellCountSum[cell[3]] + cellCount[cell[3]]] = n;
            //cout << "OK" << endl;
            ++cellCount[cell[3]];
            //cout << "OK2" << endl;
            //std::cout << r[0] << "  " << r[1] << "  " << r[2] << std::endl;
            //std::cout << cellCountSum[cell[3]] << std::endl;
            //if (n == 10) {
            //    exit(1);
            // }
        }
        //cout << "OKOKOK" << endl;

        Arrays.fill(NN, 0);

        for (int n1 = 0; n1 < number; ++n1) {
            final double[] r1 = {x[n1], y[n1], z[n1]};
            findCell(r1, cutoffInverse, numCells, cell);
            for (int k = -1; k <= 1; ++k) {
                for (int j = -1; j <= 1; ++j) {
                    for (int i = -1; i <= 1; ++i) {
                        int neighborCell = cell[3] + (k * numCells[1] + j) * numCells[0] + i;
                        if (cell[0] + i < 0)
                            neighborCell += numCells[0];
                        if (cell[0] + i >= numCells[0])
                            neighborCell -= numCells[0];
                        if (cell[1] + j < 0)
                            neighborCell += numCells[1] * numCells[0];
                        if (cell[1] + j >= numCells[1])
                            neighborCell -= numCells[1] * numCells[0];
                        if (cell[2] + k < 0)
                            neighborCell += numCells[3];
                        if (cell[2] + k >= numCells[2])
                            neighborCell -= numCells[3];

                        for (int m = 0; m < cellCount[neighborCell]; ++m) {
                            final int n2 = cellContents[cellCountSum[neighborCell] + m];
                            if (n1 < n2) {
                                double x12 = x[n2] - r1[0];
                                double y12 = y[n2] - r1[1];
                                double z12 = z[n2] - r1[2];
                                double[] newXYZ = applyMic(x12, y12, z12);
                                x12 = newXYZ[0];
                                y12 = newXYZ[1];
                                z12 = newXYZ[2];
                                final double d2 = x12 * x12 + y12 * y12 + z12 * z12;
                                if (d2 < cutoffSquare) {
                                    NL[n1 * MN + NN[n1]++] = n2;
                                    if (NN[n1] > MN) {
                                        System.out.println("Error: number of neighbors for atom " + n1 + " exceeds " + MN);
                                        System.exit(1);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }


    public void run(String fileName) {
        try (BufferedWriter br = new BufferedWriter(new FileWriter(fileName))) {
             //时间积分过程
             toOVITO("animation.xyz", 0.0);
             for (int step = 0; step < numSteps; ++step) {
                 if (neighbor_flag != 0) {
                    //如果使用了邻近列表技术，则找邻近原子, 其中已经设置了周期边界条件
                    findNeighbor();
                } else {
                    //如果什么方法都不使用，则设置周期边界条件
                    applyPbc();
                }
                //时间积分
                integrate(true, timeStep);
                //求解力
                computeLJForce();
                //时间积分
                integrate(false, timeStep);

                if (step % Ns == 0) {
                    //计算动能
                    final double kineticEnergy = computeKineticEnergy();
                    //由动能算温度
                    final double T = kineticEnergy / (1.5 * K_B * number);
                    //计算温度, 动能, 势能
                    br.write(T + " " + kineticEnergy + " " + pe + "\n");
                    System.out.println("Step " + step + " of " + numSteps + " steps");
                    if (step != 0) {
                        toOVITO("animation.xyz", step * timeStep);
                    }
                }
            }
            System.out.println(numUpdates + " neighbor list updates");
        } catch (IOException e) {
            e.printStackTrace();
        }

    }

    public void toOVITO(String fileName) {
        try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
            bw.write( number +"\n");
            bw.write("Lattice = \"");
            bw.write(boxMatrix[0][0] + " " + boxMatrix[0][1] + " " + boxMatrix[0][2] + " "
                    + boxMatrix[1][0] + " " + boxMatrix[1][1] + " " + boxMatrix[1][2] + " "
                    + boxMatrix[2][0] + " " + boxMatrix[2][1] + " " + boxMatrix[2][2] + "\" ");
            bw.write("Properties=species:S:1:pos:R:3\n");
            for(int i = 0; i < number; i++) {
                bw.write(type[i] + "  ");
                bw.write(x[i] + "  " + y[i] + "  " + z[i] + "\n");
            }
        } catch (IOException e) {
            e.printStackTrace();
        }
    }

    public void toOVITO(String fileName, double time) {
        if (Math.abs(time) < 1.0e-10) {
            try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
                bw.write(number + "\n");
                bw.write("Lattice = \"");
                bw.write(boxMatrix[0][0] + " " + boxMatrix[0][1] + " " + boxMatrix[0][2] + " "
                            +boxMatrix[1][0] + " " + boxMatrix[1][1] + " " + boxMatrix[1][2] + " "
                            +boxMatrix[2][0] + " " + boxMatrix[2][1] + " " + boxMatrix[2][2] + "\" ");
                bw.write("Properties=species:S:1:pos:R:3  ");
                bw.write("Time = 0.0\n");
                for(int i = 0; i < number; i++) {
                    bw.write(type[i] + "  ");
                    bw.write(x[i] + "  " + y[i] + "  " + z[i] + "\n");
                }
            } catch (IOException e) {
                e.printStackTrace();
            }
        } else {
            try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName, true))) {
                bw.write(number + "\n");
                bw.write("Lattice = \"");
                bw.write(boxMatrix[0][0] + " " + boxMatrix[0][1] + " " + boxMatrix[0][2] + " "
                        +boxMatrix[1][0] + " " + boxMatrix[1][1] + " " + boxMatrix[1][2] + " "
                        +boxMatrix[2][0] + " " + boxMatrix[2][1] + " " + boxMatrix[2][2] + "\" ");
                bw.write("Properties=species:S:1:pos:R:3  ");
                bw.write("Time = " + time + "\n");
                for(int i = 0; i < number; i++) {
                    bw.write(type[i] + "  ");
                    bw.write(x[i] + "  " + y[i] + "  " + z[i] + "\n");
                }
            } catch (IOException e) {
                e.printStackTrace();
            }
        }
    }

    public static void main(String[] args) {
        MolecularDynamics mdModel = new MolecularDynamics(3);
        mdModel.initializePosition("xyz.in");
        mdModel.initializeTimeIntegrationInfo("run.in");
        mdModel.initializeUniformVelocity(-1, 1);

        long startTime = System.nanoTime();
        mdModel.toOVITO("InitialConfiguration.xyz");
        String outputFile = "thermo.out";
        mdModel.run(outputFile);
        mdModel.toOVITO("FinalConfiguration.xyz");
        long endTime = System.nanoTime();
        System.out.println("Elapse Time: " + (endTime - startTime) / 1e9 + " s");

    }


    private int nDim;
    private int number;
    private double[][] boxMatrix;
    private double[][] boxMatrixInv;
    private double[] thickness;
    private double pe;
    String[] type;
    double[] mass;
    double[] x0, y0, z0;
    double[] x,  y,  z;
    double[] vx, vy, vz;
    double[] fx, fy, fz;

    int[] NN, NL;
    ArrayList<Integer>[] neighbor_list;
    private final int Ns = 100;             // 输出频率
    private final double K_B = 8.617343e-5; // 自然单位下的Boltzmann常数
    private final double TIME_UNIT_CONVERSION = 1.018051e+1; // from natural unit to fs
    private int numSteps;
    private double timeStep;
    private double temperature;

    private final int MN = 1000;    //近邻列表中允许存在的最大原子个数
    private final double cutoffNeighbor = 10.0;      //近邻列表截断半径
    private int neighbor_flag = 1;     //使用近邻列表指标器
    private int numUpdates = 0;      //记录更新近邻列表次数
}